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Large networks of coupled oscillators appear in many branches of science, so that the kinds of 
phenomena they exhibit are not only of intrinsic interest but also of very wide importance. In 
1975, Kuramoto proposed an analytically tractable model to describe these systems, which 
has since been successfully applied in many contexts and remains a subject of intensive 
research. Some related problems, however, remain unclarified for decades, such as the 
existence and properties of the oscillator glass state. Here we present a detailed analysis of a 
very general form of the Kuramoto model. In particular, we find the conditions when it can 
exhibit glassy behaviour, which represents a kind of synchronous disorder in the present case. 
Furthermore, we discover a new and intriguing phenomenon that we refer to as super- 
relaxation where the oscillators feel no interaction at all while relaxing to incoherence. Our 
findings offer the possibility of creating glassy states and observing super-relaxation in real 
systems. 
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The Kuramoto model (KM) 1 was introduced and developed 
to provide an analytically tractable description of the 
populations of coupled phase oscillators that so often 
appear in real life. It has been applied successfully in many 
fields 2 ' 3 , for example, to describe the collective behaviour of 
lasers 4 ' 5 , neurons in the brain 6 ' 7 , Josephson junction arrays 8 and 
even humans 9 . 

The widespread applications of the KM to real problems has 
ensured that the basic model, together with its variants and 
modifications 10-24 , has been studied very extensively over the 
past few decades. The practical usefulness of these studies stems 
from widespread abundance of user- defined systems that can be 
described by the KM, for example, the above-mentioned laser 
arrays 4 ' 5 and Josephson junction circuits 8 , where one can adjust 
the coupling between the oscillators quite freely. Therefore, if 
some new behaviour has been uncovered theoretically, it can 
immediately be implemented and observed in practice. 
Furthermore, for those systems where the exact configuration is 
not known precisely and cannot be changed, such as interacting 
neurons in the brain 6,7 , it is obviously useful to know what kinds 
of dynamics the different forms of KM can demonstrate, so that 
the observed behaviour can be matched with an appropriate 
existing model and thus explained. 

The KM has been found to exhibit a remarkable diversity of 
interesting phenomena and states and, by now, most of them 
have been thoroughly investigated. There is, however, one 
important exception — the so-called 'oscillator glass' state 25 . 
Thus, many systems, including networks of interacting 
spins 26,27 , dipoles 28 , and electrons 29 , have been found to 
display behaviour reminiscent of a glass structure 30 , and it has 
been suggested 25 that populations of coupled oscillators can also 
demonstrate glassy states of some kind. Studies of such states are 
expected to be very fruitful, in the same way as were the studies of 
spin glasses, which led to many new techniques and applications 
in other fields (biology, computer science, economics and so 
on 27 ). However, the provenance and properties of the oscillator 
glasses are still subject to debate 25,3 1 . 

Here we report the analysis of a very general, yet analytically 
tractable, form of KM. We derive a set of equations describing its 
steady- state behaviour and show that, for a particular class of 
distributions, the consideration may be reduced to a much 
simpler model. Thus, in the limit £— >oo, the system can 
demonstrate the same macroscopic behaviour for different 
coupling configurations, enabling one to select and study the 
simplest one. This is, in itself, an interesting theoretical result. It 
also implies that there may be cases where it is, in principle, 
impossible to infer the underlying coupling structure from the 
observed data. 

We also discuss when and how the full-time evolution of the 
system parameters can be obtained and find some interesting 
features related to this question. Among the phenomena that the 
model can exhibit, we find states with a glassy structure that can 
be studied analytically within the framework presented, and we 
discuss these states and their properties in detail. Finally, we 
describe a completely new phenomenon which we refer to as 
super- relaxation where, under certain conditions, the coupling 
between the oscillators effectively disappears during their 
relaxation to incoherence. 



Results 

Outline of the results. To make what follows clearer and easier to 
understand, we start by summarizing the main results of the 
work. First, we derive equations (18-22) describing the stationary 
states (SSs) of the system (1). Then we reduce the macroscopic 
steady- state behaviour of the system (1) (five distributed 



parameters) to a much simpler model (29) (two distributed 
parameters) for a large class of distributions (equation (23)). We 
discover new states in the coupled oscillator populations, 
which are in some sense similar to physical glasses, and we 
establish the conditions for their appearance (equations (36) and 
(37)). And finally, we discover the emergence of a super- relaxa- 
tion phenomenon, where the oscillators evolve interaction -free, 
irrespectively of the couplings between them, if the conditions 
(41), (42), (43) are fulfilled. After reviewing the background and 
terminology, each of these results is obtained and considered in 
detail below. 

The model. We consider the KM in the form 

k- N 

Oi coi -^^2% sin(0 f - Oj - ft - yj), (1) 
;=i 

where N is the number of oscillators, 0 f (£) and a>i are the fth 
oscillator's phase and natural frequency, respectively, and kflj 
(Pi + 7j) represent the coupling strengths (phase lags) between the 
fth and jth oscillators; all parameters T={(D b k h q b fi b y t } are 
drawn from a joint probability density G(T)=G(cd, k, q, y). 
The case considered earlier 10 corresponds to (1) with g,-=l, 
Pi = Vi = 0- 

The KM has not previously been treated in such a general form 
(1), though it actually includes as special cases the many KM 
modifications and extensions studied earlier. Thus, not much is 
known about the possible behaviour of the system except in those 
particular cases. To study it, we first generalize the recently 
presented framework 10 to encompass (1), and then we proceed to 
a consideration of the phenomena that it can exhibit. 

Main equations. The oscillators' collective behaviour in (1) can 
be described by two complex parameters 

z - - A Y - ™* - ^f> iT 'A (2) 

j=l j=l 

where Z is the mean field whose amplitude £ quantifies the extent 
of the agreement between the oscillators' phases 6 h while Y 
represents the weighted mean field, with amplitude W reflecting 
the agreement between 0; + y z ( + 7i for ^<0), weighted by 
\q\. With the use of the definition (2), the model (1) can be 
rewritten as 

0* = o)i - kiWsm(6i - ft - <D). (3) 

As can be seen, the dynamics of the system is governed mathe- 
matically by the weighted mean field 7, which determines the 
effective interaction between the oscillators. The ordinary mean 
field Z, on the other hand, represents a more physical macro- 
scopic variable, being the one that is most often observed in 
practice (for example, it quantifies the overall output current of a 
Josephson junction array 8 ). Note that, unlike the mean fields' 
phases O and *F individually, their difference O — *F is invariant 
under the phase shifts — (p(t), and thus represents another 

meaningful parameter. 

In the continuum limit AT— > oo, the system (1) is treated using 
the probability density function /(0,r,£), which reflects the 
probability that the oscillator has parameters T and phase 9 at 
time t. It can be further factorized as 

f(6,T,t) = p(B,t\T)G(T) (4) 

where p(6,t\T) is the conditional probability density function 
(CPDF), reflecting the probability that at time t the oscillator has 
a phase 9, given its parameters T. By definition, the CPDF should 
satisfy p_ % p(0, t \ T)d6 = 1, which leads to the continuity 
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equation 

d t p(6,t\T) + d e {[co-kWsin( 



-/j-o)]p(e,t|r)} = o, (s) 



where we have used expression (3) for 6. 

The CPDF can usually 37-40 be represented using the ansatz 
introduced by Ott and Antonsen (OA-ansatz) 41 



p(0,t\T) 



l+2Re 



/" 



a(T,t)e- w 
l-a(T,t)e- 

j e p(0,t\T)dO = a(T,t), 



1 

2n 



|a(r,oi<i 



(6) 



where we use a(T,i) = cc*(T,i) instead of the original OA- 
variable 41 a(r,£)> and the * denotes complex conjugate. 

Finally, by substituting equation (6) into equations (5) and (2), 
one obtains the full system of equations describing the dynamics 
of the system (1): 

da - ■ kW < ? -KP + Q) _gf(0+<D)) = q 5 (7) 



dt 



icoa H (a 

2 v 



We 1 



^fl(r,t)G(r)dr, (8) 



Re 1 " 



J a(T,t)G(T)dT, (9) 

where here and below, unless otherwise specified, the integration 
over dT=dcodkdqdpdy is taken over the whole domain (( — oo, 
oo ) for co, k, q and ( — n, n) for y). 

Parameter redundancy. It should be clarified that parametrizing 
the KM (1) with both p { and yj is, in a strict sense, mathematically 
redundant, since one of these phase shifts can be removed by a 
change of variables. Thus, in terms of 9\ — 9i — P iy the system 
equation (1) becomes 

^ = ^ - ^Z) 9 sin ^- - e J -V> y> = yi + Pi ( 10 ) 

with the new distribution of parameters G(T) being given as 
G(co,k,q,y) = J G(co,k,q, P,y)S(y — P — y)dpdy. Taking into 
account the relationship between the variables in equations (10) 
and (1), it is clear that the conditional distribution of the new 
phases p(9,t\co,k,q,y) is related to the original one as 
p(0,t\co,k,q,(i,y) = p(Q — ft,t\(D,k,q, ft + y). Substituting this 
into equation (2), it can be shown that the weighted mean fields 
in terms of 0 f and are equal, Y(t) — Y(t), but that the ordinary 
mean fields can, in general, change in a non-trivial fashion, 
Z(t) / Z(t). 

Nevertheless, both p t and yj can sometimes be practically 
relevant, that is, have physical meanings. For example, such a 
coupling structure can be designed artificially in the case of user- 
defined systems (for example, laser arrays 4 ' 5 ), but we expect it to 
appear in the observation-only systems as well (for example, 
interacting neurons 6,7 ). Hence, in some cases, the full model (1) 
will actually provide a more straightforward and meaningful 
description of the system, while equation (10) will represent only 
a mathematical formulation. Thus, for example, Z might be 
related to the real physical quantity,_such as the total power 
output of the laser array 4 ' 5 , while Z will then be a purely 
mathematical characteristic. To preserve generality, therefore, and 
to avoid complicating the discussion by introducing different 
variable transformations, we study the model (1) including both P 
and y. Furthermore, as will be seen below, the distributed P also 
has significant consequences in terms of the non- equilibrium 
dynamics, which cannot be straightforwardly obtained from the 
transformed system equation (10). 



It should also be noted, that the model (1) is invariant under 
(ki,Pi) -> ( - kiyfc + ri) or (q b yd -> ( - qtfi + n) or rescalings 
(ki,qj)^>(ki/r,rqi); to preserve generality, we retain this ambiguity, 
as the normalization can be fixed at any time, and there is no 
universal choice (see Methods). 

Terminology and notation. We adopt terminology similar to 
that introduced earlier 10 ' 42 . The KM form (1) is invariant under 
9' = 9 — Qf, which just changes the natural frequencies to 
co' -co — Q, so the parameter distribution becomes G'(a/, k, q, 
fl, y) — G{(d' + Q, k, q, p, y). Thus, one can consider the KM in 
different frames rotating at frequency Q with respect to some 
reference frame. For the latter, we select the frame with zero mean 
frequency < co ) = J coG(T)dT — 0 and call it the natural frame; 
the distribution G(T) is defined in this frame. We define a SS as a 
state with time-independent CPDF d t p(0,t\T) = Q, which also 
implies d t Z — d t Y—0. The state can be stationary only in a 
particular rotating frame, so it is characterized by its frame 
frequency Q and mean fields Z,Y. We refer to SSs with Q = 0 as 
natural states, while SSs with Q ^ 0 are travelling wave states. For 
later convenience we define 

G± (r) = G( ± co + Q, k, q, p, y), L(x, A) - A/ K 



P{4>,r) 



x 2 +A 2 ' 

(l-r 2 )/2n r- x (\-r 2 )e^ jln 



(1 +r) 2 +4rsm 2 ((j)/2) (C* -r)(j- 1 - e^) ' 



(11) 



where 0<r<l; the distribution P(0,r) and the way to simulate 
it are discussed in Methods. Note, that P(0,O) = l/2n and 
P(0,l) = 5(0). 

We also introduce 

I(co,k) = J e ifi G(T)dqdpdy, l£ = I( ± co + Q, k), 

J(co,k) = J qe^ + ^G(T)dqdpdy, J± = J(±co + Q,k). 

(12) 

As will be seen below, I and / represent an effective complex 
distribution of co,k to be used for the determination of Z and 7, 
respectively. 

Validity of the OA-ansatz. In what follows, we will make fre- 
quent use of the OA-equations (7-9). For the case k { — const, 
qi—ly Pi = }'i = 0 and for a very large class of frequency dis- 
tributions, their validity has been proven 39 ' 40 in the asymptotic 
limit f->oo, corresponding to system's steady- state behaviour; 
though not justified rigorously, this result seems to hold for 
distributed k, q, y as well (based on numerical simulations). 
As the next step, Pikovsky and Rosenblum 37 ' 38 derived more 
general equations and showed that the full dynamics of the 
model— for all t and almost any parameter distribution— obeys 
equations (7) -(9) only if the initial phase configuration also 
belongs to the OA-manifold 



p{9,0\T) = P{9 - <j) 0 {T),r Q (T)) 



(13) 



with any cj) 0 (T) and r 0 (T). Although not given explicitly by 
Pikovsky and Rosenblum 37 ' 38 , equation (13) can be deduced from 
their equation (3) in ref. 38, which generates (13) in the case of a 
uniform distribution of the constants of motion \jj k (for which the 
OA- description was proven to hold 37 ' 38 ); see also the work of 
Marvel et al. 43 Note that the OA-ansatz often provides a good 
approximation to the system dynamics even when equation (13) is 
not satisfied (for example, for p(0,0\T) = R(0)d(6) + (1 - R(0))/2n). 
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In addition, it might sometimes be useful to analytically 
continue a(T,t) into the complex plane over some of T, for 
example, to consider co to be complex. This is required to apply 
the conventional OA-reduction procedure 41 (see also refs 13-15) 
which, where possible, allows one to obtain a finite-dimensional 
system of equations for Y(t),Z(t). However, one should always 
have |a(r,t)|<l as, otherwise, the OA-ansatz (6) becomes ill- 
defined. Hence, a(T,t) can be considered only in the region of T 
for which 



Kr,o)| = 



J e'V(0,O|r)d0 



< 1, 



d t \a(T,t)\ 



w=1 =-Im(co) + (W/2)x 

Re[fc(e'V^-e-'V ) ')] <0 



(14) 



(15) 



at any W and (p=arg[a(r,f)] — O. The condition (14) establishes 
that \a\ < 1 is satisfied at t = 0, while (15) then guarantees that it 
holds at all other times too. Some manifestations of the issues 
related to (14) and (15) can be found in Section 3.2.0.1 of ref. 38 
and in refs 11,13-15. Note that (14) is always satisfied when there 
is no correlation of initial phases with the system parameters T: 



p(6,0\T)=P(6-*¥(0),R(0)). 



(16) 



Unless otherwise specified, we will assume that the system starts 
from such a configuration. 

Stationary states. Having reviewed the background and related 
issues, we now proceed to the derivations. We start by finding 
possible configurations into which the system (1) can settle as t-> 
oo . Although in specific cases it can converge to some inherently 
time- dependent solution, such as a standing wave 43 or an 
oscillating n- state 44 , we restrict the consideration to SSs only, 
treating them in frames where they are stationary. The effective 
parameter distribution therefore becomes G(T) -> (r), with Q 
denoting SS frame frequency. By definition and equation (6), the 
SSs satisfy d t a = 0. Using this in equation (7) and taking account 
of the OA validity condition \a\ < 1, one finds that all possible SSs 
in their own rotating frames correspond to 



fls (D = e i (° + «x. 



VMM^+m if |a,|<|jfc|w, 



• co — sign(eo) y/co 2 — k 2 W 2 



kW 



if I col > \k\W. 



(17) 



The stationary distribution of the oscillators' phases is then 
recovered by substituting equation (17) into equation (6), which 



gives 



arcsin(^) —0 + nH( — k))if \co\<\k\W, 



i ^Pfa if \co\>\k\W : 

\a> — kW sin(t/ — p — ®J| I I I I ' 



(18) 



where H( — k) denotes a Heaviside function. As can be seen, the 
oscillators with \co\ < \k\ W (in a rotating frame) are frozen around 
the positions determined by co, k and the weighted mean field 
Y— We*®, while the others are incoherent. 

It should be noted that, in addition to a s (T) given by equation 
(17), there exists one more stationary solution of equa tion (7): the 
same as equation (17), but with a minus before the V k 2 W 2 — co 2 
for |cu|<|fc|W; however, this solution is never realized in reality 
because it corresponds to the unstable position on the phase 
circle, as can be seen by recovering the respective CPDF (which 
will be the same as equation (18), but with + n for 

M<|Jfc|W). 



Self-consistency and stability conditions. While the solution 
(18) gives a qualitative picture of how the phases are distributed 
in the stationary regime, one cannot infer directly from p s (9\T) 
the stability of such a configuration or the macroscopic para- 
meters W,Q for which it can be realized. To find the latter, we 
first note that the equation (17) was derived from equation (7) 
only, without taking account of equation (8), which it should also 
satisfy. Therefore, substituting equation (17) into equation (8), 
one obtains the self- consistency conditions (SCCs) 



-\k\W 



Vk 2 W 2 -w T / fi + dco + i J coJ + dco 



--my 

-i r Vca 2 -kW[J+ -7 n -]d«l = F[Jq , W], 

J \k\W J 

(19) 

where Jq are as defined in equation (12). Taking the real and 
imaginary parts of equation (19) yields two equations from which 
the SS parameters W,Q can be determined. These can then be 
used to determine Z as well, for which, using equation (17) in 
equation (9), one gets 



Re 1 



,/(¥-<D) 



(20) 



The SS stability, on the other hand, is generally hard to study 
analytically. For the incoherent state (W^O), however, it can be 
analysed using the approach of ref. 45. Thus, performing a linear 
stability analysis of (5) above incoherence and using a self- 
consistency argument (see Methods), one can show that 
incoherence changes stability when there exists a solution x to 
the equations 



' KlmJkJ(x,k)dk + RtJ k ^±^ dcodk = 0, 
^ rcRe JkJ(x, k)dk - lmf k,{m + x ' k) dcodk = 2. 



(21) 



Interestingly, from equations (21) and (12) it follows that in the 
case G(r) = G(co, q, y)p(k) the stability of incoherence does not 
depend on a particular form of p(k), but only on the mean 
coupling strength {k}=fkp(k)dk y as noticed previously for a 
simpler model ' 15 . Note also that, if phase shifts y or ft are 
present then, with increasing coupling strength, the incoherence 
can not only lose, but also gain 46 ' 47 ^ stability at the transition 
points determined from equation (21). 

To estimate (at least approximately) the stability of the SSs with 
W>0, one can utilize the approach of ref. 10 and devise the 
empirical stability conditions (ESCs) which, for the considered 
model (1), take the form (see Methods) 



d w ReF +W(dnImF) < 1 , 
(0 w RqF - l)(d Q lmF) - (d w lmF)(d n ReF)>0. 



(22) 



where F = F[Jq , W] (see equation (19)). Despite being empirical 
and thus approximate, the ESCs (22) work well in the majority of 
cases, though not in all, for example, they can fail in the presence 
of standing waves 10 ; their performance also becomes less good 
when arg[/(co,A;)] ^0. Nevertheless, ESCs seem to be exact if 
arg[/(co,/c)] = 0 and the distribution G(T) is unimodal over co. 

Uncoupled distributions. As discussed above, for the system (1), 
the parameters of its possible SSs can be found from the SCCs 
(19) and (20), their stability can be deduced from equation (21) 
for the incoherent state, and (approximately) from equation (22) 
for other SSs, while the associated phase distributions are given by 
equation (18). The corresponding expressions, however, are 
generally quite complicated, but it turns out that they simplify 
greatly if we consider the distribution of q, y to be uncorrected 
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with a>, k. In what follows, we will therefore assume that 

G(co,k,qJ,y) = g(co,k)h(qJ,y), g£ = g(±co + Q,k), 

(23) 

so the effective (cl>, k) -distributions (12) become 

I(co,k) =\I\^"g(<o,k), J(a>,k) =\J\<»'g(a>,k), 



|/| e >4>, = 



J e'Phiq J,y)dqdpdy, 
\J\j*' = y qe i(fi + y) h(qj,y)dqdpdy. 



(24) 



The SCCs (19) then simplify to 

F w (W,0) = Wcas<l>j/\J\, F Q (W,Q) = - Wsin^/ (25) 
where 

dfc 

-|fc|w 



F W (W,0) = — g+VkW-ctfdco, 

J KW J _\ k \ w 
*W«)-/^{/«*b + *»- (26) 



-£; 2 W 2 da; 



\k\W 



) 



(28) 



are fully analogous to F# 5 q in ref. 10; the expression (20) for £ and 
O — reduces to 

R = \I\W/\f\, <b-V = <j> J -4> I , (27) 
the incoherence stability conditions (21) become 

\/^(Q,fc)dfc = 2 ^, 

and the ESCs (22) can be simplified using 

F[Jn,W]=\J\j*j(F w + iFci)- 

System reduction. From equations (25-28) it is evident that, for 
uncoupled distributions (23), all the macroscopic properties 
of the SSs are completely characterized by g(a>, k), |/| and 4>j (see 
equation (24)), irrespectively of the particular form of h(q, y), 
while |J|, 0/ serve merely to specify Z (27). Instead of (1), 
therefore, one can consider the system 



COi - 



W\<sr ■ m 



(29) 



with the same distribution of a> h k { defined by g(a>, k). Then each 
SS of (29) corresponds to an SS of (1), and their parameters are 
related as 



Q = Q SK , W=|/|#sk, iH'I^SK, 



(30) 



where the subscript SK denotes the states of the model (29); the 
stability of the corresponding states will also be the same, as 
follows from equation (28) for incoherence, and from the ESCs 
(22) and numerical evidence for other SSs. 

Hence, for any distribution of q, y obeying (23), one can 
reduce the consideration of the steady- state behaviour of (1) to 
the much simpler Sakaguchi-Kuramoto model (29). This elegant 
and unexpected result is illustrated in Fig. 1. Interestingly, 
something similar was noted earlier 21 , but for a significantly less 
general model than (1). The reduction (30), however, relates only 
to the macroscopic properties of SSs (f-> oo), whereas the full 
evolutions Z(t),Y(t) and the microscopic properties of the 
resultant states cannot be obtained in this way. Note that 



0.4 
Hf 0.2 




O Original system 

0 Reduced system (rescaled) 



10 
K 



15 



20 



Figure 1 | Numerical confirmation of system reduction. The SS 

parameters W (green), R (red) and Q (black), as calculated for the original 
system (1) (circles) and obtained from the reduced system (29), (30) 
(diamonds), are shown in dependence on coupling K; solid lines correspond 
to theoretical predictions based on (25-28). The distribution for the original 
system was:g(o;,/c) - 3{k - K)[oj 2 + e"^] ~\h{qj, y) - e -(q-V 2 /o.02 
[P(fi - tc/8, 0.8) + P(p + tc/2, 0.9)]P(y - tt/3, 0.7), but the same picture 
appears for any h(q, ft, y) with identical 4>j, \J\, \l\ (24). The simulations 
were performed for 500 s, and the values presented are averages over the 
last 100 s. 



equation (30) implies W< \J\,R<\I\, that is, the distribution of q, 
y imposes an upper bound on the SSs' mean field strengths, 
which they cannot exceed however strong the coupling is. 

Definition of glassy states. Ordinary glass 30 is a state of matter 
which, in contrast to the liquid state, is solid; but, in contrast to 
the crystalline state, has a structure without any long-range 
translational order. States that are in some sense similar were 
reported in several other systems 26-29 , but whether or not there 
can exist glassy analogues in networks of coupled oscillators — 
oscillator glasses — has remained unclear 25 ' 31-36 . 

However, a very basic question is how to define an oscillator 
glass, that is, what properties the system should possess in order 
to qualify for this title. Obviously, one should draw an analogy 
with known glassy states in order to illuminate this issue. Thus, 
the defining feature of the structural glasses 30 — frozen disorder- 
can be specified as: (i) the particles, for example atoms, lack any 
long-range translational order, similarly to the liquid state; (ii) the 
particles are frozen, that is, do not move with respect to each 
other, similarly to the solid state. Later, states with frozen spins 
lacking long-range orientational order were discovered in spin 
systems, and by analogy were called spin glasses 26 ' 27 . The 
canonical spin glasses, in addition to properties (i) and (ii), are 
characterized by: (iii) a significantly redundant ground state, so 
that the same macroscopic behaviour can be realized by many 
microscopic spin configurations, not related to each other by any 
simple symmetry transformation; (iv) frustration, that is, no spin 
configuration can simultaneously satisfy all energy bounds; 
(v) slow non-exponential relaxation of the order parameter 
(magnetization for spin systems), as well as other exotic 
dynamical features. 

Relating the phases Q { of the oscillators to the positions of 
the particles in space, or the spin orientations, the analogies 
of the above properties for populations of oscillators will be 
(all for f-> oo ) 



Uniform distribution of phases at any time : 
p(0,t) = /p(0,t|r)G(T)dr=l/27t; 

Existence of a (frozen) group of oscillators 
not moving with respect to each other : 

Q = WuhT)E \(e^- 6 ^)\>0, 



(31) 



(32) 



where (...) denotes time average; 
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Nonuniqueness of lim^oo — Oj(t)] 
for the frozen oscillators, which might , . 

therefore end up in different places, ^ ' 

depending on the initial conditions; 

Frustration, that is, (arguably) nonexistence 

of stable phase configurations with all (34) 

0 f (f) = 0, even if all a>i = 0; 

Slow nonexponential relaxation of R(t): , x 

Mm^R^ = oo, W > 0. ^ ] 

Generally, it is quite a subtle issue to decide what is a 'true' 
oscillator glass, that is, should it satisfy all the criteria (31-35), or 
only some of them? Thus, the state reported by Daido 25 has 
properties (31), (34) and (35), but not (32) and (33), since the 
oscillators there, although adjusting their mean phase velocities to 
some extent, perform a diffusive motion and thus, are never 
frozen (Q = 0). This state, therefore, bears more analogy to spin 
glasses than to structural glasses. 

Here, we utilize a different approach, considering state to be 
glassy if it satisfies the first two criteria (31), (32). To avoid 
possible terminological issues, we refer to such states as quasi- 
glassy. They are thus denned as the states with a uniform 
distribution of phases 6 t (31), indistinguishable from incoherence 
but where, in contrast to the latter, some macroscopic portion of 
oscillators are frequency-locked (32), that is, have equal phase 
velocities fy. 

It should be noted that the condition for the absence of phase- 
locking is sometimes taken as £ = 0. However, although being 
implied by equation (31), it is not as rigorous as the latter. Thus, 
there might be many configurations for which R = 0 but where 
the oscillators actually adjust their phases. Examples include two 
phase-locked populations of the same size, which are in the anti- 
phase with each other. States with R — 0 satisfying condition (32), 
but not condition (31), will be called spurious glassy, according 
to ref. 31 where they were first discovered. All SSs with R>0, 
including the usual synchronized states, n- states 19,20 and 
travelling waves, will be classified as coherent. 



As a simplified picture, the distinctions between the states can 
be understood in terms of a group of people doing cyclical 
exercises, each with their own tempo and other parameters. The 
incoherent state is when everyone proceeds independently; the 
coherent state is when they all move synchronously, at any given 
time having similar poses; spurious glassy is when they exercise at 
the same tempo but remain pairwise in the opposite poses, and 
quasi-glassy is when they adjust their tempos, but always remain 
in the independent random poses, thus representing a kind of 
synchronous disorder. 

Realization of glassy states. Considering model (1) with the 
uncoupled distribution (23), it can be shown (see Methods) that if 
the marginal distribution of /3 is uniform and |/| > 0, then all SSs 
except incoherence satisfy equations (31) and (32). Therefore, the 
quasi-glassy states appear when 

J h(qj,y)dqdy = l/2n, (36) 

W>0^\J\= | J qe Ky + ^h(qJ,y)dqdl]dy\>0. (37) 

Representing h(q, 7)=/z 1 (/?)/z 2 (^, y\P), equation (36) becomes 
hi(P) = 1/271. Thus, to satisfy also condition (37), the distribution 
of q, y should be specifically correlated with The simplest 
examples are h(q, y\ft) — S(q — q 0 )P(y + /?,r) and h(q, y\ft) — 
L(q - ^ocosjft, A)P(y - y 0 > r). 

For spurious glassy states, the condition (36) is not satisfied, 
but R = 0, which in the present case, is equivalent to 
\I\=\fe^hi(P)dp\ = 0, as can be seen from equation (27). There 
are many possible /zi(/?) satisfying this condition, for example, 
h 1 (P) = P(P-p 0 ,r)+P(ft-p 0 -n,r). 

Examples of the four types of states occurring in the model (1) 
are shown in Fig. 2, presented in the natural frame where the 
population does not move as a whole (<&>> =0). In the (6,cd)- 
plane, the only difference between the quasi-glassy state and 
incoherence is that, in the former, phases within the glassy cluster 
(\co\ < \k\W) are frozen around random angles (static disorder), 
as seen in (0, /?) -plane; this is in contrast to their asynchronous 




Figure 2 | Different kinds of states. The states that can appear in the model (1), illustrated by snapshots of positions in the 6, co (upper panel) and 6, ft 
(bottom panel) planes of (the same) 1,000 randomly selected oscillators (out of N = 25,600). Tiny vertical arrows show the oscillators movements during 
0.25 s (for a dynamic version of the figure see Supplementary Movie 1). Red-dashed lines show boundaries between the clusters and incoherent 
populations. In all cases, g(co, k) = L(co, l)d(k-K), h 2 (q,y\P) = <5(q -1)<5(y + A and (a) K = l fc,(j5) = 1/2n; (b) K=3, fc,(j9) = 1/2n; (c) K=3, M#) = 
(1/2)[P(j8-7c/2,0.8) + P(j8 + 7c/2,0.8)]; (d) K=3, h^P) = P(p, 0.8). The pictures will be the same for any h 2 (q, y\fi) satisfying \J\=\ </>j = 0 in each case. 
All states have Q = 0 and are presented in the natural frame (<a>> =0); depending on the distribution G(r) of parameters in (1), similar states 
may also appear in a rotating frames (Q#0), corresponding to quasi-glassy, spurious glassy or coherent travelling waves. 
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movement (dynamic disorder), observed for incoherence. Draw- 
ing an analogy between the oscillators' phases and the particles' 
positions in space, the incoherent state of the system (1) can be 
related to a liquid or gaseous state of matter; coherent to the 
crystalline state; spurious glassy to a crystal with a specific 
structure; and quasi-glassy to a glass. 

As discussed previously, the model (1) can be reduced to a 
form without /3 (10) by a change of variables. Thus, any quasi- 
glassy or spurious glassy state in terms of 0 Z - ^1) can be 
transformed to a coherent state in terms of 0 Z - = 0{ — $ { . 
Contrariwise, glassy states of the system (1) can always be viewed 
as coherent states but in a specifically disordered coordinate 
system, being therefore redundant in a strict mathematical sense. 
In practice, however, the coordinates are chosen based on their 
physical meaning, so that these states are realizable when the 
actual interaction between the oscillators contains phase shifts 
which can be achieved, for example, in a real systems with 
programmable coupling 48 ' 49 . A similar situation occurs for the 
Mattis spin glass 50 : it can be transformed to a ferromagnetic state 
by change of variables 26,27 . Nonetheless, the Mattis model found 
a number of applications 51-53 and proved to be a useful starting 
point in studies of spin glasses. Note that, combining the results 
of the present work with those of refs 22,32,33,54, one can modify 
the model (1) to observe states exhibiting most of the glassy 
properties (31-35) (see Methods). 

Relaxation dynamics. Up to now, we have concentrated mainly 
on the behaviour of the system (1) in the asymptotic limit t-> oo , 
restricting our consideration to a set of its possible SSs. The 
relaxation to these states is determined by the parameter dis- 
tribution and can take various forms. It is hard to analyse in 
general, but for particular G(T) the full-time evolutions Y(t) ,Z(t) 
can be obtained analytically by using the OA-reduction proce- 
dure 41 , which can be extended to incorporate the distributed 
phase shifts. Thus, consider the distribution 

G(T) = L(co, A)S(k - K)S(q - l)S(fi)P(y - 0 O , r 0 ) (38) 

In this case, a(T,t) can be analytically continued inside the unit 
circle of z y = e iy , as the condition (15) is satisfied there. Hence, in 
equations (8) and (9), one can integrate over y by changing 
the integration to the unit circle of z y and taking the residue at 
the pole z y = r 0 e^° of P(y — 4> 0 ,r 0 ) (see equation (11)). The 
integration over co is performed in the usual way 41 , that is, by 
taking the residue at cd — /A. Then, from equations (8) and (9) 
one obtains Z(t) = a(/A, e iy = r 0 £*, t), Y(t) = r 0 £*Z(f). 
Substituting this into equation (6) and solving the resultant 
equations, one gets finally 



W(t) 



l/l 



2A| 
K, I 



^i + [(i-p/|/w(o)) 2 -i] e -«-^| 



4>(f) = Q(t) 



Kj tan</> 



W 2 (t) 



+ 1 



Kj = JC|/|COS0; 



R(t)=\I\W(t)/\J\, <D(f)-<F(f) 



(39) 

where we have represented all in terms of J, / (see equation (24)), 
which in the present case are |J|e z ^ 7 = 1, |/|^ = r 0 e^ 0 . 
Obviously, a similar procedure to that outlined here can be 
applied if the distribution of y has a few poles inside the 
unit circle. 

Now consider the KM (1) with distributed 

G(T) = L(fl), A)d(k - K)8{q - l)P(p - <j> 0 , r 0 )5{y), (40) 
in which case, one has Z(t) = Y(t) by definition (2). By the change 



of variables 9i = 6i — /? z , the present case can be reduced to the 
already- studied KM withjiistributed y (see equation (38)). Then, 
because Z(t) — Y(t) = Y(t) (see discussion of equation (10)), 
the macroscopic parameters should still evolve according to 
equation (39), although one should now use \ l\e 1 ^ 1 =|/|e 2 ^ = 
r 0 e^° in the latter. The same result may be obtained by applying 
the OA-reduction, with the integrals (8) and (9) over being 
evaluated by taking residues inside the unit circle of e^, similarly 
to what was done in the previous case. 

However, the behaviour predicted by equation (39), and the 
actual behaviour of Z(t),Y(t), agree only for distributed y, but not 
for distributed as demonstrated in Fig. 3. The reasons for this 
are, first, that one cannot continue a(T,t) inside the unit circle of 
as (15) is not satisfied there, which makes OA-reduction 
impossible when /? is distributed. Second, the evolution of 
Y(t),Z(t) for (40) cannot be obtained from the system of 
0j = 6j — at least if starting from the standard initial 
conditions (16), as assumed here. This is because such 
transformation introduces a correlation of the initial phases with 
system parameters. For example, if all 0/(0) = 0, then one will 
have 0i(0) = -ft = -% so that 2(r,0) = e~ iy and the OA- 
reduction cannot be applied because now (14) is not satisfied. 
Therefore, although the cases (38) and (40) can be transformed 
into each other by changes of variables, the time evolution for 
distributed /? appears to be much more complex than for 
distributed y, and cannot be obtained in a simple form. 

Super-relaxation. Astonishingly, for a class of distributions G(T), 
and a very large family of initial configurations, the oscillators do 
not feel any interaction at all while relaxing to incoherence. This 
behaviour occurs when the following three conditions are satisfied: 



/ 



qe ty G(co, k, q, y)dkdqd/3dy — 0, Vcl>; 



Initial conditions p(0, t = 0 1 T) are not 
correlated with q or y, either directly or 
indirectly (through the other parameters) 

Incoherence is the only stable state 
of the system. 



(41) 



(42) 



(43) 



Thus, conditions (41) and (42) imply W(0) = fe i6 p(6,0\T)qe iy G(T) 
d6dT = 0. Based on numerical evidence, the weighted mean field 
then stays at zero during the whole evolution W(t) = 0, leading to 
an effective disappearance of interaction between the oscillators, as 
follows from equation (3). As a result, the oscillators evolve freely 
(0; = (Di) and so the relaxation (both its rate and form) depends 



DC 0.5 - 




Time (s) 

Figure 3 | Theoretical and numerical system behavior. The main plot and 
inset display time evolutions of R(0 and W(0, respectively. Dashed lines 
show the behaviour predicted by (39), while solid lines correspond 
to the results of numerical simulations: blue— distributed y (38), red- 
distributed p (40). In (38) and (40), we used (j) 0 = n/6, r 0 = 0.5. 
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g(cD,k)=L(coS) P (k), h(q,M=S(q-mMV2n 



p(k)=S(k) 
. p(k)=d(k-V 
■ p(k)=u(k-1,1) 



— 1 



g((o,k)~e-^d(k-^, h(q,p,y)~p(q,y\P)[P(P-2n/3,0.8)+P(P+n/3,0.8)] 




p(q,y\P)=S(q-cosP)P(y-n/6,0.8) 
. . > p(q,y\P)=u(q^S)P(y+P-7i/6,0.8) 
p(qf,y| J 6)=L(qf-cos J 6,0.1)u(y-7r/6,7i/4) 



g((D,k)=u{(D^)8(k-^,h{q,M=S(q-^m^n 




PW)=W) 
p(P)=V2n 

p(p)=0.3P(p-n/6,0.5)+0.7P(p-n/2,0.8) 



Time (s) 



10 



Figure 4 | Examples of super-relaxation. Interaction-independent relaxation of the order parameter to incoherence for different parameter distributions 
(all satisfying (41)), which are indicated in the figure; u(x, b) denotes uniform distribution of x in [ — b, b], while L(x, A) and P(x, r) are defined in (11). 



Incoherence (0<K< 2) — K=5 — K=-\0 — K=20 




0 1 2 3 4 5 6 



Time (s) 

Figure 5 | Two-stage relaxation to quasi-glassy states. Time evolutions of 
R(0 (solid) and W(t) (dashed) from the initial conditions R(0) = 1 to quasi- 
glassy states (g(co, k) = L(co, W(k-K), h(q, p, y) =S(q -1)5(0 + y)/2n) 
for different constant couplings /c, = /C. The inset shows the results on a 
logarithmic ordinate scale. 



only on the marginal distribution of co: 

R(t) = K(0)| J ^G(r)dr|. (44) 

This phenomenon, which we refer to as super-relaxation, is 
illustrated in Fig. 4. Note that for uncoupled distributions, 
equation (23), the condition (41) reduces to Jqe iy h(q y y) 
dqd/Jdy = 0. 

The claim that W(t) = 0 if all the conditions (41), (42) and (43) 
are satisfied can be proven rigorously for particular distributions, 
for example, G(T) = G(co, k, q, /3 + y)/2n (see Methods), though 
we were unable to prove it in general; based on simulations, 
however, it seems to be satisfied in most (if not all) cases. Thus, it 
is clear that the system always has a solution with W(t) — 0 if 
conditions (41) and (42) are fulfilled: in this case, to have 
W(t) >0, the distribution p(6, t\T) should acquire a particular 
dependence on the system parameters; but evolving according to 
equation (5) with W(t) = 0, it can become dependent only on co, 
which in turn cannot increase W(t) above zero if condition (41) 
holds. The stability of this solution is rather hard to prove, but it 
appears to be stable if the incoherence is the only stable state 
(condition (43)). Nevertheless, even if W(t) is not exactly zero, 
below some relatively high threshold its value seems to have a 

8 



negligible effect on the dynamics of R(t), as will be seen in Fig. 5 
below. 

Super-relaxation can appear only if the final state is incoherence 
(W^O), while relaxation to SSs with W>0 will generally be 
coupling dependent, even if conditions (41) and (42) are satisfied. 
In the latter case, relaxation occurs in two stages as shown in Fig. 5 
for the example of quasi-glassy states. The phases first begin to 
disorder in the same way as for incoherence, but are slowly 
entrained while passing their equilibrium positions. When the field 
of the entrained oscillators, characterized by W, becomes strong 
enough, they begin to force the unentrained ones to take their 
positions, so the relaxation switches to a faster, coupling- dependent 
regime. This switch occurs sooner for stronger coupling. 

Discussion 

We have generalized our earlier approach 10 to make it applicable 
to the more general KM (1) so that, using equations (19-22), one 
can immediately obtain a macroscopic characteristics of possible 
SSs. The generalized KM (1) encompasses a variety of KM 
variants studied earlier 10 ' 13-15 ' 19 ' 21 ' 24 , allowing one readily to 
reproduce and extend many of the previous results. Remarkably, 
the steady-state behaviour of (1) with any distribution h(q y p y y) 
(see equation (23)) can be obtained from the simple Sakaguchi- 
Kuramoto model (29), (30). It should be noted, however, that 
(17-22) describe only SSs, being inapplicable to inherently 
non- stationary solutions such as standing waves 55 or oscillating 
7i-states 44 . Note also that most expressions can straightforwardly 
be extended to the case G(T) G(T y R y W) (see Methods); for 
examples when this might be relevant see refs 56-59. 

Most interestingly, we have found, that the model (1) can 
exhibit exotic behaviour, such as glassy states and super- 
relaxation, thereby opening new horizons for KM-related 
investigations, both theoretical and practical. These discoveries 
have a far-reaching implications. For example, it should now be 
possible to create, observe and study glassy behaviour in real 
systems of coupled oscillators, where a variety of novel 
phenomena may be anticipated. As one possible application, if 
some physical quantity can be associated with the weighted mean 
field W (equation (2)) in laser arrays 5 described by the KM, it 
might be possible to construct a laser exhibiting zero intensity 
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CR = 0) but for which the other effects are nonvanishing (W>0); 
similar considerations apply to a wide range of different KM 
applications. Furthermore, the phenomenon of super-relaxation 
might be used to design systems whose dynamics remains highly 
stable in the face of different perturbations and parameter 
changes. 

Methods 

Numerical simulations. Except where otherwise specified, all simulations were 
performed by integrating the full system equation (1) using a 6th order Runge- 
Kutta method with a time step of 0.01 s. The number of oscillators used was 
N = 25,600 for Fig. 2,N= 10 5 for Fig. 1 and N= 10 6 for all other cases. Additional 
details are provided in figure captions. 

Parameter normalization. Since the model (1) is invariant under (k b /?,•)->( — k b 
ft + n) or (q b y { ) -> ( — q b y { + n) or rescalings (k b q { ) -» (fc,-/r, rqD> there exist some 
ambiguity in the parameter definitions. To remove this ambiguity, one can fix the 
normalization of q { and specify a rule for choosing the sign of k b q b but different 
choices might be preferred in different situations. For example, the normalization 
J2<li = 1 is inapplicable when q t = cos(y/) with y,- being uniform; = 1> on tne 

other hand, fails for a Lorenzian distribution of q b The most universal choice seems 
to be f\fqeW + v ) G(T)dqdpdy\dcodk = 1, which additionally assures W<1. How- 
ever, when q { = 1 and ft are distributed (so Y = Z), this normalization will lead to a 
rescaling of q { (so Y=£Z), which might be inconvenient. Therefore, because the 
normalization can be fixed at any time, we retain the associated ambiguity in order 
to preserve generality. Note that W (equation (2)) but not R and not \k\W, changes 
under the {k,q) -rescalings and thus can be higher than unity. 

Phase distribution. The distribution P(cf), r), equation (11), represents a Poisson 
kernel for the unit disc and is very common for the KM. For example, it 
can be shown that the OA-ansatz (equation (6)) can be rewritten 43 as p(9, 
t\T) = P(0 — arg[a(I\ t)], \a(T, t)\). Because P(</>, r) has one pole inside and one 
outside the unit circle of z$=i^, it is very convenient in relation to analytic 
derivations, for example, one has f n _ n e^P(</> — 4> 0 , r)d</> = re 1 ^ 0 . To simulate this 
distribution, one sets 0 f = 2 arctanfj^tan^p, — 1/2))], where p, e [0,1] are 
uniformly distributed random numbers. 

Derivation of the incoherence stability conditions. To perform a linear stability 
analysis of incoherence, we examine the maximum growth rate for small pertur- 
bation £rj(0, t\T) added to the incoherent solution p(9, t\T) = 1/271. Following 
ref. 45, this perturbation is expressed as 



unified form for all co, k as 



rj(9,t\T) = c(T)e~ l0 e At + c* (I>V f + rj ± (0, 1 1 r), 



(45) 



where the Fourier expansion of r] ±(9, t\T) contains only terms e in ® with \n\ > 1 (so 
that j(e ±i6 )rj ± (9, t\T)d0 = 0). Substituting p(9,t\T) = l/2n + eri(9,t\T) and (45) 
into (5), and collecting the terms proportional to e~ l ®, to the first order in £ one 
obtains 



Ec(T)(k-ico)e M 



Y(t). 



(46) 



Considering the part ~ e l of equation (5), one will get the complex conjugate of 
equation (46), while for terms with \n\ > 1 the growth rate of the associated 

perturbations Q t [sr] ±(9,t\T)] can be shown to be 0(e 2 ), so that they are of no 
interest for linear analysis. 

From (46) it follows that c(T) = {ke^Y{t)e~ Xt )/(2s(k - ico)). Substituting this 
into Y(t) = f e i9 p(9,t\T)qe iy G(T)d9dT = se h J 'q^ 'c(T)G(T)dT and dividing both 
sides by Y(t)/2, one self- consistently obtains 



/' 



kqe ^ + ^G{Y)dY _ 
k — ico 



J 



kj(m, k)dmdk 



(47) 



where we have used the definition of J(co,k) (see equation (12)). The incoherence 
changes stability when Re[k] crosses zero. Hence, denoting k = k T -\- ik v taking the 
limit l r — >0 in equation (47) and using lim^ r ^ 0 [^r — i{ m ~ ^0] ~ 1 = n ^{°° ~ h) 
+ i[oo — ki]~ l , one obtains (21) with x = k^. 

Derivation of the ESCs. The ESCs (22) cannot be derived rigorously but, rather, 
are based on empirical assumptions. Thus, we first assume that the perturbations to 
the SS obey 



SY = A 



f 



a s (T,Y + SY)G+ +xl (T)dT-(Y + 5Y) 



(48) 



where dY= (SW+ iWS<&)e and 3Q, are the deviations of the SS parameters from 
their stationary values Y and Q, respectively; A is a constant, Gq is defined in (11) 
and a s (Y,Y) is the OA stationary solution (17). The latter can be represented in 



«.(r,ir) = 



jp V\k 2 \Yf -a) 2 \exp[-i arg[fc2|7 j 2 " ffl2] ] + to 



kY* 



(49) 



Although equation (48) represents a purely empirical assumption, some of the 
motivation behind it is that for incoherence it gives the correct transition points 
(equation (21)), as can be checked by performing a linear stability analysis of 
equation (48) with 7=0. 

For 7>0, based on trial and error we set A>0, 30. = W 2 (5<D in equation (48). 
Then using equation (49) and retaining in (48) only the terms of the first order in 
3W, <5®, one obtains 



3W 



W~ 



ReF-1^ 
l {d w lmF) 



W 2 (d a ReF) 
W(dalmF) 



( 5W 



(50) 



where we have expressed all in terms of F = F[Jq , W] (equation 19). For A>0, 
the system (50) is stable only if the trace and determinant of the corresponding 
matrix are lower and higher than zero, respectively, which gives (22). 

Glassy conditions for the model considered. For the system equation (1) with 
the uncoupled distribution, equation (23), it is easy to see that the conditions 
(31,32) for observing quasi-glassy states can be rewritten as (36,37). Thus, the 
phases of the oscillators with |co z | < \ki\W in a stationary regime are frozen at 
9j = ^ + arcsin(cL),/|/c ; | W) + nH( — k { ) (all in the rotating frame), as follows from 
equation (18). Therefore, for states with W>0, an absence of phase locking 
between the oscillators (condition (31)) is equivalent to a uniform distribution of /? z - 
(equation (36)). Next, W>0 (37) is the necessary and sufficient condition for (32): 
necessary because, otherwise, all oscillators have = co f = 0, as implied by 
equation (3); and sufficient because it establishes frequency locking of the oscil- 
lators with \a>i\ < \ki\W (since their phases are constant). Finally, as is clear from 
equation (25), for the uncoupled distributions (23), SSs with W>0 can appear only 
if |/| >0, as additionally indicated in equation (37). 

Glassy states exhibiting nearly all glassy properties. One can modify the model 
(1) to observe glassy states with most of the desired properties (31-35), and not 
only the first two of them. Thus, frustration (34) can be introduced, for example, by 
considering the van Hemmen type of interactions 32 ' 33 , while the high multiplicity 
of states (33) can be achieved by using the higher-order coupling 22 ' 54 . As an 
example, it can easily be checked numerically that a state satisfying conditions 
(31-34) appears in the system 



Y^Kii + qtkj) sin[3(0 f - dj - ft - jj)], (51) 



where k { and q { are independent random variables taking values + 1 with equal 
probability, ft = — y { are uniformly distributed and the frequencies CO; are drawn 
from g (co) ~ [1 + co 2 ] ~ l . Regarding the last glassy property (35), it might be 
possible to satisfy it with an appropriate choice of the parameter distribution, 
which, to a large extent, determines the relaxation of R(t). However, most 
interestingly, it turns out that for systems of the type (51), there exist many 
configurations with R(t^> oo) >0, and the latter is proportional to the initial #(0). 
Thus, in many cases the phases do not disorder completely, so the basic criterion 
(31) becomes violated if starting from initial conditions other than incoherence. 

Super-relaxation claim lV(t) = 0. In some cases, one can rigorously prove that 
W(t) remains zero at all times if the conditions (41), (42) and (43) are fulfilled. 
Thus, consider the distribution 



G(r) = G(co,k,qJ + y)/2n. 



(52) 



Changing the variables to 9[ = 0,- — ft, one obtains the system (10) with 
G(r) = G(ct), k, q,y), characterized by the same weighted mean field W(t) = W(t) 
(see discussion of (10)). Under such a change, however, one subtracts from each 9 t 
a uniform random variable ft, so that any initial conditions for 9j satisfying (42) 
will be mapped to a uniform initial condition for 9[\ moreover, since G{T) depends 
only on the sum y + ft these initial 9[ will not be correlated witlijy,- = y { + ft, as 
would be otherwise. Therefore, one has p(0,-, 0 1 Y) = l/2n =>• a(T, 0) = 0. Then, 
taking into account that W(0) = W(0) = 0, it follows from equations (7) and (8) 
that a(f , t) = 0 => W(t) = W(t) = 0. Finally, because the incoherent state for 9 { 
is stable (condition (43)), it will obviously be stable in terms of 0; as well, so that the 
perturbations to the weighted mean field will decay, ensuring that it stays at zero. 

Generalization to G(r)->G(r f R f W). In some cases 56-59 , the parameter 
distribution might depend on R and/or W, for example, the coupling can be 
influenced by the mean field strength G(Y, R, W) ~ 3(k — K(R)). One can easily 
extend all expressions describing the SSs to such cases. Thus, proceeding in the 
usual way, the stationary OA solution (17) and the associated phase distribution 
(18) can be shown to have the same form. Next, the SCCs equations (19) and (20), 
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from which the SS parameters are determined, will also be the same except 

I(co, k) -> I(co, k, R,W) = / e ip G(T, R, W)dqdfidy, 

J r (53) 
J(co,k)^J(co,k,R,W) = / qe^ p + 7) G(T,R,W)dqdpdy, 

Note, however, that if the parameter distribution depends on R (d R G(T,R,W)^0), 
and q or y are distributed (so that W=£R), then equations (19) and (20) become 
coupled and should be solved at the same time to get all the R, W,Q,0 — x ¥. Such 
situation occurs because, in this case, the effective interaction between the 
oscillators, although explicitly determined by W (3), has additional implicit 
dependence on R; as a result, the dynamics of the model becomes more 
sophisticated, and its analysis appears to be an interesting topic for future research. 

Assuming that lima_> 0) w->o R n W m d R l d^G(T 1 R, W) = 0, Vn, m > 0, n + m > 0, 
the conditions for incoherence stability can be derived in the usual way and take 
the same form (21), except that /(o>, k)^f(co, k, 0, 0) (53). The ESCs (22), on the 
other hand, cannot so easily be generalized to field-dependent distributions, at least 
if d R G(T y R,W)^0 and W^R; otherwise, that is, when G(T)^G(T, W), they 
preserve the original form of equation (22), but with a modified / given by equation 
(53) (though it is not clear how well they work in this case). Having obtained 
the general equations, their simplification for the uncoupled distributions 
G(r, R, W) =g(co, ky R y W)h(q, jff, y, R y W) to analogues of equations (24-28) 
is straightforward. 
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